function pso = MinPSO(fn,range,amount,speed,trgt,maxStep)
  %% 基于粒子群算法的优化算法(Particle Swarm Optimization)
  %
  % fn 目标函数
  % range 参数范围 [x1_Min,x1_Max;x2_Min,x2_Max;...]
  % amount 粒子数
  % speed 速度权重 [惯性权重,自身最优权重,全局最优权重]
  % trgt 目标值
  % maxStep 最大步数
  %
  %% 计算公式:
  %
  % 全局最优:X' = [x'1,x'2,x'3....]
  % 自身最优:X'' = [x''1,x''2,x''3....]
  % 自身当前:X = [x1,x2,x3....]
  % 前速度:V' = [v'1,v'2,v'3....]
  % ----------------------------------------------
  %
  % 现速度:V = V'*speed[1] + (X''-X)*speed[2] + (X'-X)*speed[3]
  %         = speed * [V';(X''-X);(X'-X)]
  % 现当前:X = X+V
  %
  % 计算适应度 fn(X)
  %
  %% 数据结构
  %
  % X_Best<x,fx> 全局最优解
  % XList = struct[]('X','V','Best','BestFx')
  %
  LastBest = [];
  [XList,best] = MinPSOInit(fn,range,amount,speed);

  for i=1:maxStep
    if (i>1) && (LastBest - best{2} <= trgt)
      fprintf('%d\n',i);
      break;
    end
    LastBest = best{2};
    [XList,best] = MinPSOProcess(fn,XList,amount,speed,best);
  end
  pso = best;


end

function [outList,best] = MinPSOInit(fn,range,amount,speed)
  %% 初期化粒子群
  %
  % fn 目标函数
  % range 参数范围 [x1_Min,x1_Max;x2_Min,x2_Max;...]
  % amount 粒子数
  % speed 速度权重 [惯性权重,自身最优权重,全局最优权重]
  %
  % outList 粒子群 struct('X'List,'V'List,'Best'List,'BestFx'List,'AllBest'{x,Fx})
  %

  % dimensionality(维度)
  dim = size(range,1);
  outX = [];
  outV = [];
  % 随机生成 粒子群 初始位置和速度
  for i=1:dim
    outX = [outX,rand(amount,1)*(range(i,2)-range(i,1))+range(i,1)];
    outV = [outV,(rand(amount,1)*2-1)*sum(speed)];
  end

  fxList = cellfun(fn,mat2cell(outX,ones(amount,1)));
  [bestFx,bestFxIndex] = min(fxList);

  best = {outX(bestFxIndex,:),bestFx};
  outList = struct('X',outX,'V',outV,'Best',outX,'BestFx',fxList);

end  % function

function [outList,best] = MinPSOProcess(fn,inList,amount,speed,best)
  %% 计算粒子
  %

  % 现速度:V = V'*speed[1] + (X''-X)*speed[2] + (X'-X)*speed[3]
  %         = speed * [V';(X''-X);(X'-X)]

  speedFn = @(v,x,bestX)speed*[v;bestX-x;best{1}-x];
  vList = arrayfun(speedFn,inList.V,inList.X,inList.Best);

  % 现当前:X = X+V
  xList = inList.X+vList;

  % 计算适应度 fn(X)
  fxList = cellfun(fn,mat2cell(xList,ones(amount,1)));

  updateList = inList.BestFx > fxList;
  inList.BestFx(updateList) = fxList(updateList);
  inList.Best(updateList,:) = xList(updateList,:);

  inList.X = xList;
  inList.V = vList;

  [bestFx,bestFxIndex] = min(inList.BestFx);
  best = {inList.Best(bestFxIndex,:),bestFx};

  outList = inList;
end
